	#clean up
	rm(list=ls())
	
	#set seed
	set.seed(715130425)
	
	#load libraries
	library(akima)
	library(geoR)
	library(scatterplot3d)  
	library(ggmap)
	library(lattice)
	library(foreign)
	library(maps)
	library(spdep)
	#install.packages('/Volumes/rgdal/rgdal_0.9-1.tgz',repos=NULL)
	library(rgdal)
	#install.packages('/Users/jamie/Downloads/prevR_2.9.tgz',repos=NULL)
	library(prevR)
	library(car)
	library(xtable)
	library(maptools)
	
	#options
	options(scipen=8)
	
	#Set to working directory. Adjust the command below to reflect YOUR working directory.
	#setwd('/Volumes/MONOGAN/gillAreal/ccesCode/')
	#setwd('/Users/jamie/Documents/gillAreal/ccesCode/')
	#setwd('/Users/jamie/Desktop/nhgis0003_shape/nhgis0003_shapefile_us_tract_1970/')
	
	#world map
	a<-map("world",regions="USA",resolution=0)
	dev.off()
	b<-cbind(a$x,a$y)
	#write.csv(b,"map50.csv",col.names=F,row.names=F)
	#b<-read.csv("map50.csv")
	b[b[,1]>0&!is.na(b[,1]),1]<- -180-(180-b[b[,1]>0&!is.na(b[,1]),1])
	plot(x=b[,1],y=b[,2],pch=".")
	points(x=-117.1611,y=32.7157,pch=3,col='red') #San Diego, CA
	points(x=-147.7164,y=64.8378,pch=3,col='red') #Fairbanks, AK
	points(x=-124.73,y=48.16,pch=3,col='red') #Cape Alava, WA, westernmost point on mainland
	points(x=-122.4781,y=48.7491,pch=3,col='red') #Bellingham, WA
	points(x=-122.6765,y=45.5231,pch=3,col='red') #Portland, OR
	points(x=-134.4197,y=58.3019,pch=3,col='red') #Juneau, AK
	points(x=-149.9003,y=61.2181,pch=3,col='red') #Anchorage, AK
	#abline(v=-130,col='blue') #easternmost point in Alaska in the western hemisphere
	#abline(v=-141,col='blue') #Alaska-Canada border
	points(x=-118.2437,y=34.0522,pch=3,col='red') #Los Angeles, CA
	points(x=-122.4194,y=37.7749,pch=3,col='red') #San Francisco, CA
	points(x=-157.8583,y=21.3069,pch=3,col='red') #Honolulu, HI
	points(x=-155.681111,y=18.9111,pch=3,col='red') #Ka Lae, southernmost point in Hawaii
	points(x=-154.811,y=19.516,pch=3,col='red') #Capa Kumukahi Light, easternmost point in Hawaii
	points(x=-124.4095,y=40.4401,pch=3,col='red') #westernmost point in California
	points(x=-121,y=34.0522,pch=3,col='red') #Due west spot of Los Angeles
	abline(h=50,col='blue')
	abline(h=23,col='blue')
	b<-cbind(b,as.numeric(b[,2]>50))
	b<-cbind(b,as.numeric(b[,2]<23))
	points(x=b[b[,3]==1,1],y=b[b[,3]==1,2],pch=".",col='blue')
	points(x=b[b[,4]==1,1],y=b[b[,4]==1,2],pch=".",col='green')
	
	#reproject
	b.2<-na.omit(b)
	b.3<-SpatialPoints(b.2[,1:2],proj4string=CRS("+proj=longlat"))
	b.4<-spTransform(b.3, CRS("+proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=37.5 +lon_0=-96 +x_0=0 +y_0=0 +ellps=GRS80 +datum=NAD83 +units=km +no_defs"))
	b.5<-cbind(b.4@coords,b.2[,3:4])
	
	plot(b.5,pch=".")

	#figure out new movements
	movements<-rbind(c(-141,61.2181),
	c(-124.73,32.7157),
	c(-154.811,21.3069),
	c(-124.73,45.5231))
	movements.2<-SpatialPoints(movements,proj4string=CRS("+proj=longlat"))
	movements.3<-spTransform(movements.2, CRS("+proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=37.5 +lon_0=-96 +x_0=0 +y_0=0 +ellps=GRS80 +datum=NAD83 +units=km +no_defs"))
	movements
	movements.3@coords
	
	#move the states
	c<-b.5
	c[c[,3]==1,2]<- c[c[,3]==1,2]-(3456+348.6)+500 #Move Alaska
	c[c[,3]==1,1]<- c[c[,3]==1,1]-(-2287+ 2337)
	
	c[c[,4]==1,2]<- c[c[,4]==1,2]-(217+348.6)+450  #Move Hawaii
	c[c[,4]==1,1]<- c[c[,4]==1,1]-(-6009 + 2337)+75
	
#	png("AKHIreproject.png",width=960,height=480)
#	postscript("AKHIreproject.eps",width=8,height=4,paper="special")
	plot(x=c[,1],y=c[,2],pch=".",axes=F,xlab="",ylab="",asp=1,lwd=2,xlim=c(-5482,2255),ylim=c(-1331,1633))
	points(x=c[c[,3]==1,1],y=c[c[,3]==1,2],pch=".",col='red')
	points(x=c[c[,4]==1,1],y=c[c[,4]==1,2],pch=".",col='blue')
#	dev.off()
	
	
###Plot a ZCTA vs. ZIP code overlay###
###Load ZIP code locations from Tom Tom data###
#Points accessed on September 2, 2015 from http://www.arcgis.com/home/item.html?id=1eeaf4bb41314febb990e2e96f7178df
#Polygons accessed on September 2, 2015 from https://www.arcgis.com/home/item.html?id=8d2012a2016e484dafaac0451f9aea24
zip.polys.0<-readShapePoly(fn="Zipcode_Polys/Zipcodes_Polys.shp",proj4string=CRS("+proj=longlat +datum=NAD83"))
#zip.polys.1<-zip.polys.0[!zip.polys.0$STATE%in%c("AK","GU","HI","PR","VI"),]
#zip.polys.1<-zip.polys.0[zip.polys.0$STATE%in%c("GA"),]
zip.polys.1<-zip.polys.0[zip.polys.0$ZIP_CODE%in%c(30601),]
plot(zip.polys.1)
zip.30601<-as.matrix(zip.polys.1@polygons[[1]]@Polygons[[1]]@coords)

	###ZCTA map###
	#Downladed July 2, 2016 from https://www.census.gov/geo/maps-data/data/cbf/cbf_zcta.html
	zcta.polys.0<-readShapePoly(fn="cb_2014_us_zcta510_500k/cb_2014_us_zcta510_500k.shp",proj4string=CRS("+proj=longlat +datum=NAD83"))
	#zcta.polys.1<-zcta.polys.0[!zcta.polys.0$STATE%in%c("AK","GU","HI","PR","VI"),]
	zcta.polys.1<-zcta.polys.0[zcta.polys.0$ZCTA5CE10%in%c(30601),]
	plot(zcta.polys.1,lty=2)
	zcta.30601<-as.matrix(zcta.polys.1@polygons[[1]]@Polygons[[1]]@coords)
	
	#draw the overlay
	#postscript("zip30601.eps",width=6,height=6)
	#pdf("zip30601.pdf",width=6,height=6)
	plot(x=zcta.30601[,1],y=zcta.30601[,2],type='n',axes=F,xlab="",ylab="")
	polygon(x=zcta.30601[,1],y=zcta.30601[,2],lty=2,border="red",lwd=2)
	polygon(x=zip.30601[,1],y=zip.30601[,2],lty=1,border="blue",lwd=2)
	#points(x=-83.308837,y=34.011828)#my house
	legend(x=-83.308,y=34.05,col=c("red","blue"),lty=c(2,1),lwd=2,legend=c("ZCTA","ZIP"))
	#dev.off()

###ANOTHER EXAMPLE OF A ZCTZ VS. ZIP OVERLAP FROM STANFORD, CA###
zip.polys.1<-zip.polys.0[zip.polys.0$ZIP_CODE%in%c(94305),]
plot(zip.polys.1)
zip.94305<-as.matrix(zip.polys.1@polygons[[1]]@Polygons[[1]]@coords)

	zcta.polys.1<-zcta.polys.0[zcta.polys.0$ZCTA5CE10%in%c(94305),]
	plot(zcta.polys.1,lty=2)
	zcta.94305<-as.matrix(zcta.polys.1@polygons[[1]]@Polygons[[1]]@coords)
	
	#draw the overlay
	#postscript("zip94305.eps",width=6,height=6)
	#pdf("zip94305.pdf",width=6,height=6)
	plot(x=zcta.94305[,1],y=zcta.94305[,2],type='n',axes=F,xlab="",ylab="")
	polygon(x=zcta.94305[,1],y=zcta.94305[,2],lty=2,border="red",lwd=2)
	polygon(x=zip.94305[,1],y=zip.94305[,2],lty=1,border="blue",lwd=2)
	legend(x=-122.195,y=37.445,col=c("red","blue"),lty=c(2,1),lwd=2,legend=c("ZCTA","ZIP"))
	#dev.off()
	
	

